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ABSTRACT 

Gamma-ray emission from pulsars has long been modeled using a vacuum dipole field. This ap- 
proximation ignores changes in the field structure caused by the magnetospheric plasma and strong 
plasma currents. We present the first results of gamma-ray pulsar light curve modeling using the 
more realistic field taken from three-dimensional force-free magnetospheric simulations. Having the 
geometry of the field, we apply several prescriptions for the location of the emission zone, comparing 
the light curves to observations. We find that when the emission region is chosen according to the 
conventional slot-gap (or two-pole caustic) prescription, the model fails to produce double-peak pulse 
profiles, mainly because the size of the polar cap in force-free magnetosphere is larger than the vac- 
uum field polar cap. This suppresses caustic formation in the inner magnetosphere. The conventional 
outer-gap model is capable of producing only one peak under general conditions, because a large frac- 
tion of open field lines does not cross the null charge surface. We propose a novel "separatrix layer" 
model, where the high-energy emission originates from a thin layer on the open field lines just inside 
of the separatrix that bounds the open flux tube. The emission from this layer generates two strong 
caustics on the sky map due to the effect we term "Sky Map Stagnation" (SMS). It is related to the 
fact that the force-free field asymptotically approaches the field of a rotating split monopole, and the 
photons emitted on such field lines in the outer magnetosphere arrive to the observer in phase. The 
double-peak light curve is a natural consequence of SMS. We show that most features of the currently 
available gamma-ray pulsar light curves can be reasonably well reproduced and explained with the 
separatrix layer model using the force-free field. Association of the emission region with the current 
sheet will guide more detailed future studies of the magnetospheric acceleration physics. 
Subject headings: MHD — pulsars: general — gamma-rays: theory — stars: magnetic fields 



1. INTRODUCTION 

Gamma-ray pulsars are of great scientific interest in 
high-energy astrophysics. The high energy emission from 
these objects provides valuable information about pul- 
sar magnetospheric structure, particle acceleration mech- 
anisms, and plasma physics in strong magnetic fields. 
Observations by the Energetic Gamma Ray Experiment 
Telescope {EGRET) onboard Compton Gamma Ray 
Observatory (CG RO) confirmed 7 gamma-ray pulsars 
( Thompson] 2004 ) . The light curves of these pulsars are 
characterized by widely separated double-peak features, 
with the first peak typically lagging the radio pulse by a 
small fraction of rotation period. Recently, A GILE mis- 
sion has detected several new gamma-ray pulsars and 
refined the detection of some previously known pulsars 



gamma-ray light curves suggests emission coming from 
the outer magnetosphere. Gamma rays are believed to be 
produced by energetic particles accelerated in the "gap" 
regions in the magnetosphere. The particles emit via cur- 
vature, synchrotron and inverse Compton (IC) radiation. 
Various theoretical models differ in the assumed location 
of the emission zones (i.e., the gaps). Conventional mod- 
els for pulsar gamma-ray emission include the polar cap 



mode l {Harding et al.||1978l iDaugherty & Harding||1982 
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19961, the slot gap model (S G, or two-pole caustic mode. 
TFUfor short; |Arons fc Scharlemann||1979HArons|[l983 



2008[ |Pellizzoni et al.||2009a]b| . More 
i year after launch, the Large Area Tele- 



Muslimov fc Harding||2003[ |2004| |Dyks~fe Ruda kj|2003 
Dyks et al .|2004j), the outer gap model(OG, |Cheng et a ' 
1986a b; Romani & Yadigaroglu|1995 Yadigaroglu|199" 



Cheng et al. 2000 ) , as well as t he inner annular gap mode 



■ 



( Halpern et al 
impressively, one 

scope (LAT) onboard Fermi Gamma-ray Space Telescope 
(Fermi) has discovered more than 40 new gamma-ray 
pulsars |Abdo et al.|2008| |2009a|c|flg|h| ). The sensitivity 
and timing capability ot Fermi- LAI have also provided 
the mos t precise light curves an d phase resolved pulsar 



spectra flAbdo et ai.]|2009b|d|e 



These new discoveries 
have greatly expanded the sample of gamma-ray pulsars 
and underscored the importance of understanding the 
origin of double-peak profiles. At the same time, the 
diversity of light curves from latest gamma-ray pulsar 
sample calls for improvements in theoretical modeling. 
The widely separated double peak feature in the 

xbai@astro.princeton.edu, anatoly@astro.princeton.edu 



(1AG, Qiao et al.|20U4| |2007[ 

Most calculations from these models approximate the 
pulsar magnetosphere by a rotating magnetic dipole field. 
In realit y, the magnetosphere is kn own to be filled with 
plasma (Goldreich & Julian 19691. This plasma is es- 
sentially TofcfriTe^fTTTT^^siynrg pE + j x B/c = 0, 
where p and j are the charge and current densities. The 
structure of the FF magnetosphere differs substantially 
from the rotating dipole field du e to the poloidal current 
and the returning current sheet (|Contopoulos et al.|1999 
Gruzinov|20"05l |Komissarov|2006| |Mckinney|iWOti| |Timo- 
khin||2006p . Recently, FF magnetospheric field structure 
in 3D has be come available from the time-dependent FF 



simulations (|Spitkovsky 2006 Kalapotharakos & Con- 
topoulos|2009 p. The FF magnetosphere satisfies E ■ B 



everywhere, hence no particle acceleration or emission 
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can formally occur. Real pulsars must lie somewhere be- 
tween the vacuum, which has E ■ B ^ everywhere, and 
FF cases. The fact that the power of gamma-ray emis- 
sion is a small, though non- negligible (< 20%), fraction 
of the total spin down power suggests that FF should 
be a reasonable approximation to the overall field struc- 
ture with the exception of relatively small regions where 
acceleration takes place. 

In the companion paper ( Bai fc Spitkovsky|2010 here- 
after BS10), we demonstrated that the modeling of 
gamma-ray pulsar light curves using the vacuum dipole 
field is subject to large uncertainties. This is because 
the peaks in the vacuum light curve are caused by the 
fortuitous overlap of radiation from different regions of 
the magnetosphere (caustics). This overlap is sensitive 
to the geometry of the magnetic field due to two effects. 
First, the shape of the field controls how the aberration 
of light and light travel delay add together to cause the 
formation of caustics. We showed that a seemingly small 
change of treating the retarded vacuum field in the in- 
stantaneously corotating frame, rather than in the lab- 
oratory frame where it should be defined, can lead to 
significant reduction in the strength of caustics in the 
two-pole caustic/slot gap model and to the moderate di- 
lution of the outer gap peaks. Second, the behavior of the 
field near the light cylinder determines the shape of the 
polar cap on the star, and, hence, indirectly controls the 
shape of the emission zone, even if it is located in the 
inner magnetosphere. We considered the sensitivity of 
the caustic formation to the modification of shape of the 
vacuum polar cap. We compared the polar cap obtained 
by tracing the vacuum field lines to a simple circular po- 
lar cap. This change significantly affected the caustics of 
the two pole caustic model. The outer gap model was 
less sensitive to this change, but, since the emission zone 
for the outer gap is in the outer magnetosphere, the field 
shape there is likely to be strongly modified by the inclu- 
sion of plasma effects. Given these uncertainties in the 
vacuum field modeling, it is highly desirable to revisit 
the existing models using the FF field configuration and 
compare the results with the latest observations. 

In this paper, we use th e force-free field fro m 3D time- 
dependent simulations by Spitkovsky (2006) (hereafter, 
S06), and present the calculation of pulsar gamma-ray 
light curves using various theoretical models of emission. 
In §2, we provide the general analysis of FF field struc- 
ture, addressing the location of the current sheet and 
its association with magnetic field lines. In §3, we de- 
scribe our numerical method and theoretical models for 
calculating the gamma ray light curves. §4-§6 are dedi- 
cated to comparison between different models. We show 
in §4 that the conventional two-pole caustic model fails 
to produce double-peak light curves using the FF field. 
We propose a novel "separatrix layer" (SL) model in §5, 
in which the emission zone is fixed in a layer just inside 
the separatrix associated with the strong current sheet. 
This layer can be accurately traced by open field lines 
originating in an annulus just inside the polar cap riirj^] 

1 In the original version of this manuscript, we termed this model 
as "annular gap" model. The name "separatrix layer" more accu- 
rately describes the location of the emission zo ne and avoids th e 
confusion with the "inner annular gap" model (Qia o et al.|2004| . 
Moreover, although the present paper is based only on geometry, 
the association of the emission layer with the magnetospheric cur- 



This model differs from the conventional TPC model in 
that the emission zone is not concentrated on the last 
open field lines, but rather interior to them, and it spans 
a larger range of heights up to and beyond the light cylin- 
der. It differs from the conventional OG model in that 
the emission zone consists of all field lines in a given 
flux tube, and not just those that cross the null charge 
surface. The double-peak light curve is a generic prop- 
erty of the FF field structure, where caustics of emission 
form in the outer magnetosphere when the field lines ap- 
proach the split-monopolar geometry. We call this "the 
stagnation effect" on the sky map. In §6, we show that 
the conventional OG model applied to FF field has diffi- 
culty in producing double-peak light curves unless a care- 
fully chosen inclination angle and viewing geometry is 
used. We discuss various implications of our results to 
the pulsar acceleration mechanisms and compare our re- 
sults with recent observations in §7. The main results 
are summarized in §8. 

2. STRUCTURE OF THE FORCE-FREE MAGNETOSPHERE 

Force-free pulsar simulations evolve a rotating magne- 
tized conducting sphere immersed in a massless infinitely 
conducting fluid. The simulations are three-dimensional 
(3D) and start with an inclined dipolar field attached 
to a sphere (S06). The resolution of our simulations 
is 80 cells per light cylinder (LC) radius, Rlc — c /^i 
and the radius of the star in simulations is 15 cells 
(i?^ m ~ 0.19i?Lc)- We ran a series of simulations vary- 
ing the inclination angle between pulsar's magnetic and 
rotation axes from a — 0° to a — 90° with an inter- 
val of 15°. The FF fields are extracted after evolving 
the FF equations for 1.2 stellar rotations. Evolving for 
longer times makes little difference to t he fiel d str ucture 



( Kalapotharakos fc Contopoulos||2009[ ). In [2.1 we il- 
lustrate the charge density and current structure of the 
FF fi eld, and present the properties of FF polar cap in 
§2.2| The 3 D cu rrent sheet structure of the FF field is 
discussed in 92.31 



2.1. Charge and Current in the Force-Free Field 

While it is not possible to determine individually the 
density and velocity of positively and negatively charged 
particles solely from the FF model, the FF current and 
charge density contain rich physical information from 
which we can gain useful insights. 

In ideal force-free electrodynamics, by assuming that 
electromag netic field patte rn corotates with the pulsar, 
we obtain ( |Gruzinov||2006"l ) 



E = —/So x B , 
V x [B + O x O x B)] = XB 



(1) 
(2) 



where fto = O x r/c is the normalized corotation veloc- 
ity, and A is a scalar function that is conserved along the 
magnetic field lines, because B ■ VA = 0. In the case of 
the aligned rotator, pa rameter A is the analog of A'(^) 
in the pulsar equation (Michel 1 973a|b ), where A is pro- 
portional to the poloidai current, and the derivative is 
with respect to the poloidai magnetic flux \I>. Therefore, 

rent sheet (see j ]7.2| implies that the acceleration may not neces- 
sarily involve a gap" at all. 
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A oc A' characterizes the current per magnetic flux. More 
specifically, positive A represents current along the mag- 
netic field line, while negative A represents current flow 
opposite to the direction of the magnetic field. A formal 
derivation for the physical meaning of A can be found in 
Appendix |A"[ where we show that AB is the field-aligned 
current density in the corotating frame. We calculate 
A by taking the dot product of B/B 2 on both sides of 
equation (pi). 

In Figure [TJ we show the color plots of A in the £1 — /i 
plane (the plane containing the rotation axis and the 
magnetic moment vector) for a range of inclination an- 
gles. We also show vectors of the magnetic field projected 
onto this plane. The strong current layer manifests itself 
as a local enhancement of A. First, consider the axisym- 
metric pulsar magnetosphere (a = 0). The field stru cture 
closely resembles the solution by |Timokhi"n| ( 2006 ) with 
x = R y/Rlc ^ 1 an d well app roaches the (JKF so- 
lution^] (Contopoulos et al. 1999). In agreement with 
these solutions, the current along the open field lines 
in our simulation is predominantly of one sign in each 
hemisphere; however, near the edge of the open zone 
there is a distributed return current region bounded by 
a strong return current sheet. Instead of a delta-function 
current sheet expected in theory, the simulation current 
sheet is several grid cells thick, giving A a finite maxi- 
mum/minimum there. The open field lines increasingly 
resemble the rotating split monopole solution beyond the 
LC and become radial in the poloidal section. The left 
panel of Figure [TJ illustrates these features. 

Next, consider the general oblique rotators. The tran- 
sition of current structure is smooth as a increases from 
0° to 90°. In Figure [TJ we see that at small inclination 
angle, a — 30° , the pattern of A is similar to the axisym- 
metric case, and the current sheet inside the LC appears 
slightly wider and weaker. As the inclination angle in- 
creases, there are also current sheet-like features inside 
the LC, but the distribution of current becomes very 
asymmetric, and more current is returned to the polar 
cap through the distributed flow rather than the current 
sheet. The current flow in the polar cap is predomi- 
nantly of one sign (e.g., ingoing) for the aligned rotator, 
bounded by a strong current sheet of opposite sign (e.g., 
outgoing), and is symmetric with respect to the equato- 
rial plane. In contrast, when a = 90°, the current flow 
has different signs in the northern and southern halves of 
the polar cap (this can be better seen in Figur e B) w hich 
shows the current through the polar cap, see j ]2?2| , and 
the current distribution is anti-symmetric with respect 
to the equatorial plane. On the periphery of the polar 
cap there are thin current layers, as shown in red (up- 
per) and blue (lower) oval structures in the right panel 
of Figure [TJ The current in these current layers forms 
loops connecting the two poles through the closed zone. 

2 In our simulation, the numerical resistivity is very low when 
the equatorial current sheet is aligned with the grid. Therefore, for 
inclination angles < 10° the evolution of the Y-point to the light 
cylinder takes longer than 2 turns to approach the expected result 
(S06). Larger inclination angles are well converged by 1.2 turns. 
We have checked this convergence by comparing results between 1.2 
turns and 2 turns (where we manage to avoid reflection by using 
a very large computational domain). No significant difference is 
found for a > 15°. In Figures 1-4, however, we do show the field 
structure after 2 turns, to allow the Y-point in the a = 0° case to 
get closer to the LC. 



Integrating over the polar cap, we find that the amount 
of current flowing to the other pole is about 20% of the 
current flowing on open field lines. Such closed current 
loops for orthogonal rotators are qualitatively consisten t 
with predictions of the model by Gurevich et al. (1993), 
although we find the amount of current shunting to the 
other pole to be smaller in our simulations. The total 
current on the open field lines (integrated by magnitude) 
is just 20% smaller for the orthogonal rotator than for 
the aligned rotator, suggesting that the current density 
for the orthogonal rotator exceeds the simple expectation 
of the speed of light times local Goldreich- Julian density 
on the polar cap. This current cannot be provided by 
charge-separated flow alone and requires abundant pair 
formation. The details of the current adjustmen t in such 
polar caps are still uncertain ( Lyubarsky|2009[ ) . 

A thin current sheet outside the LC exists for all in- 
clination angles; its structure asymptotically approaches 
the rotatin g split monopole solution outlined by |Bogov^| 
alov (1999). If the current sheet outside the LC is con- 
nected to the star, it must be connected through a cur- 
rent sheet inside the LC, because the current flow in the 
CF can n ot c ross magnetic field lines, as inferred from 
equation ( A6 ). In AppendixjB] we show that the amount 
of current m the current sheet outside the LC that is 
connected to the star monotonically decreases with incli- 
nation angle. This is related to the degradation of the 
current sheet inside the LC and the thickening of the 
Y-region for oblique rotators shown in Figure [TJ These 
features are unlikely to be caused by numericalresistiv- 
ity, and we have tested that the thickness of the Y-region 
and the strong current layers inside the LC is not sensi- 
tive to numerical resolution. For the orthogonal rotator, 
the current sheet outside the LC is totally disconnected 
from the star (see AppendixfB]). Therefore, the oval struc- 
tures in the right panel of Figure [TJ are most likely the 
current loops that connect the two poles, rather than cur- 
rent sheets connecting to the outer magnetosphere. This 
strongly contrasts with the case of the aligned rotator, 
where the current sheet inside the LC unambiguously ex- 
ists and connects to the equatorial current sheet outside 
the LC. While there are definitely strong currents flowing 
on the periphery of the open field lines for all other in- 
clinations seen from Figure [TJ determining whether their 
thickness is fini te or infinitesimal would require further 
study (see {2.2 for a speculative discussion). Through- 
out this paper, we will refer to these current sheet-like 
features inside the LC as "strong current layers" . 

The charge density necessary to provide corotation 
of the magnetosphere is referred to as Goldreic h-Julian 
(GJ) charge density (Goldreich & Julian 1969). It can 
be written as 



P : 



V • E _ U B /3 • (V x B) 



47T 



2ttc 



4-7T 



(3) 



If we assume complet e charge separation, as in |Gol-| 
dreich & Julian (1969), then the second term is of the 
order p/3qj and we arrive at the classical result with 
p cx Jl • B. Whether a real pulsar magnetosphere has 
complete charge separation is not clear, but the eclipse in 
the double pulsar system PSR J0737-3039 indicates that 
the particle density in the inner magnetosphe re is much 
higher than typical Goldreich-Julian density (Rafikov & 
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Fig. 1. — Color plots of the parameter A of the force-free field in the f! — fi plane for magnetic inclination angles a = 0° , 30°, 60° and 90° 
respectively. Arrows show the direction of the projected magnetic field in this plane. All plots have units in light cylinder radius Rlc> an d 
have the same color scale. 




Fig. 2. — Color plots of the charge density p in the force-free field in the Q — fi plane for magnetic inclination angles a = 0°,30°,60° 
and 90°. We've multiplied p by r 2 to improve the contrast. Arrows indicate the direction of the projected magnetic field in this plane. All 
plots have units in light cylinder radius Rlc> an d have the same color scale. 




Fig. 3. — Color plots of the parameter g [see eq. Q] of the force-free field in the Q — B plane for magnetic inclination angle a = 0°, 30°, 60° 
and 90°. We've multiplied g by r 2 to improve the contrast. Arrows show the direction of the projected magnetic field in this plane. All 
plots have units in light cylinder radius Rlc> an d have the same color scale. 
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2005| |Lyutikov fc Thompson||2005[ ). In Figure 



Goldreich 

51 we plot the charge density from the FF simulation 
We find that the bulk distribution of FF charge density 
is approximately proportional to £l-B. However, there is 
a significant enhancement of charge density in the strong 
current layer and in the current sheet beyond the LC, 
which is due to the contribution of the current in the 
second term of equation (|3| . 

The distribution of charge in the strong current layer 
inside the LC has another notable feature. In Figure [2j 
the charge distribution is symmetric between the north 
and south magnetic hemispheres for the aligned rota- 
tor. Dense negative charges more or less uniformly fill 
the current sheet inside the Y-point, while the current 
sheet outs ide the Y-point is positively charged (Timo- 
khin|20 06). However, this symmetry is broken when a is 
different from zero (e.g., at a — 30°). The strong current 
layer inside the LC then contains both signs of charge, 
with the northern and southern parts having opposite 
sign at fixed azimuth. The current sheet outside the LC 
also contains charges of both signs, although one sign 
of charge is prevailing (in red color). This t r end is also 
observed by Contopoulos & Kalapotharakos ( 2009 1 . 

It is also useful to visualize the b'b' magnetosphere in 
terms of space- like and time- like regions of 4-current. A 
4-current is space- like if j — p 2 > 0; otherwise, it is 
time-like. Note that this definition is Lorentz invari- 
ant. Regions with space- like current demand counter- 
streaming of different signs of charge in the frame where 
p = 0. T his could lead to plasma instabili ties and dis 



sipation (Lyubarskii 1996 Gruzinov 
sheet with space-like current could t 
produce particle acceleration and gamma -rays (iGruzinov 
2007a[ |Lyubarsky|2008[ |Spitk ovsky|2008D . In th e asymp 



2007b [ ). A current 
len be e xpected to 



totic split-monopole regime ot [5ogovalov| ( |1999[ ) , one ex- 
pects the current sheet outside the LC to be space-like. 
This is because both magnetic and electric fields change 
sign across the current sheet, and since \B\ > \E\ outside 
the current sheet, MHD jump condition demands J > o~c 
in the current sheet, where J and a are surface current 
and surface charge densitie^] In Figure [3j we plot the 
Lorentz invariant quantity 



Q = ± 



p 2 c 2 _ j2| 



(4) 



where we take plus sign for time-like current, and minus 
sign for space-like current; thus, counter-streaming is re- 
quired in the regions which are negative. We find that 
in both aligned and oblique rotators the current sheet 
outside the LC is in general space-like (blue). Although 
there are regions in red around the current sheet, the 
overall surface integrated g is negative, as expected. The 
strong current layer inside the LC appears to be pre- 
dominantly time-like. The magnetosphere in the aligned 
rotator is dominated by time-like regions, and only in re- 
gions near the null charge surface does it show very small 
negative values of g. As the inclination angle increases, 
the space-like region becomes much broader, occupying 
a substantial fraction of the open field-line regions. We 
will discuss this further in SjTJ 

2.2. The Force-free Polar Cap 

3 We thank Yuri Lyubarsky for pointing this out to us. 



The polar cap is the region on the neutron star (NS) 
surface where open magnetic field lines originate. The 
rim of the polar cap is set by the locus of last open field 
lines (LOFLs). We adopt the conventional definition of 
LOFLs, where we treat a field line as open if it crosses 
the LC; otherwise, it is regarded as closed. This is an 
oversimplification, as in FF field the closed zone does 
not have to extend to the LC everywhere; some closed 
field lines may also extend through the LC and close 
through the current sheet. However, this simple defi- 
nition facilitates comparison with calculations using the 
vacuum fieldP] In Figure |4j we plot the shape of the 
FF polar cap on the sphere of radius 0.25Rlc for incli- 
nation angles a = 0°, 30°, 60°, 90°. The surface of the 
sphere (which is larger than the NS for ease of visual- 
ization) is colored with the parameter A, indicating the 
sense of the current flow. For comparison, we also plot 
with dash-dotted lines the polar caps ob tained from th e 
retarded vacuum magnetic dipole field (Deutsch 1955). 
We note that tracing vacuum field lines beyond the LC 
always returns them back to the star, because vacuum 
solution does not have the current sheet. Nevertheless, 
as is conventional in the literature, we always call vac- 
uum field lines that go beyond the LC "open" , although 
they are still formally closed. 

From Figure [4] we see that the polar cap of the FF field 
encloses the current flowing into and out of the pulsar. 
The polar current switches from flowing into the NS at 
a = 0° , to equal halves of oppositely-directed current at 
a = 90°. The peak of the current density slightly lags 
the zero phase (i.e., the phase of the magnetic pole) due 
to pulsar rotation. The evolution of the strong current 
layer as a increases can also be clearly seen from Figure 
|4j For aligned rotator, the rim of the polar cap corre- 
sponds to the footprint of the current sheet (i.e., the red 
enhancement on the left panel), as expected. As a in- 
creases, the strong current thickens on one side (lower 
side in Figure EJ), and weakens on the other side. The 
thickened part may no longer count as a current sheet, 
as it gradually occupies half of the polar cap to become 
the main contributor of the polar current at a = 90°. 
The weakened side gradually shifts to the outside of the 
polar cap (at a — 90°), and forms the current loop in the 
closed field lines as seen in the rightmost panel of Figure 
[l] Following the discussion in the previous subsection, we 
speculate that the weakened part of the strong current 
layer actually corresponds to the current sheet inside the 
LC. The current sheet covers a circle on the polar cap for 
the aligned rotator. For oblique rotators, the region cov- 
ered by the current sheet on the polar cap reduces to an 
arc covering the upper part of the polar cap in Figure [4j 
The extent of the arc-like region shrinks with inclination 
angle a, consistent with the reduction of star-connecting 
current in the current sheet outside the LC. In the mean 
time, some current in the arc-like region on the polar 

4 Moreover, it is probably better not to distinguish between 
open and closed field lines, but to distinguish between "active" 
and "dead" field lines. Dead field lines are closed within the 
Y-point, beyond which the equatorial current sheet is launched. 
Other field lines, including open field lines and field lines that are 
closed through the current sheet, are considered active. Therefore, 
one may consider "last active field lines" as a more meaningful 
concept in the FF magnetosphere. Our definition of last open field 
lines as active field lines is equivalent to setting the radius of the 
Y-point to be at the light cylinder. 
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Fig. 4. — The shape of the polar caps in force-free field (solid lines) and the Deutsch field (dash-dotted lines) for magnetic inclination 
angles a = 0°, 30°, 60° and 90° (left to right). Overlayed are the color plots of the current to flux ratio A on the pulsar surface. To better 
visualize A, the size of the pulsar is set to be 0.25-R^o Also, the color scale of A is compressed compared to Figure^ 



cap becomes connected to the other pole, forming cur- 
rent loops. This is associated with the thickening of the 
Y-region discussed in the previous subsection. The loop 
current is zero for the aligned rotator, and increases to 
about 20% of the total polar current for the orthogonal 
rotator. 

We note that the FF polar caps are smooth and ap- 
proximately circular. On the other hand, polar caps ob- 
tained from oblique vacuum d ipole rotators disp lay notch 
structures and sharp jumps (Dyks et al. 2004). As was 
shown in BS10, gamma-ray light curves are sensitive to 
the shape of the polar caps, as this changes the forma- 
tion of caustics in the emission. The FF polar cap is also 
larger than its vacuum counterpart, due to the larger 
open magnetic flux. This is consistent with the larger 
spin down power for force-free rotator compared to equiv- 
alent vacuum rotato r (S06). The same effect is seen in 
the aligned rotator ( Contopoulos et idT| |l999 ). As a re- 
sult, the LOFLs predicted from the vacuum dipole field 
should lie within the LOFLs, or in the open volume, of 
the FF field. This effect also has significant consequences 
to the appearance of the FF sky map. 

The magnetic field lines that originate near the po- 
lar cap from NS surface can be p arameterized using 



the open volume coordinate system (Yadigaroglu 1997 
Cheng et alpOOOl |Dyks et al.||2004[ ). This system has 
two coordinates, (r ov , <p m ), recording the footprint on 
the NS surface, where <p m is the magnetic azimuth, and 
r ov is the magnetic colatitude normalized by the mag- 
netic colatitude of the LOFLs. Therefore, the LOFLs 
correspond to r ov = 1, and open field lines have r ov < 1. 

2.3. Current Sheet-Field Line Association 

Current sheets are generic features of FF pulsar mag- 
netosphere. To better visualize the current sheet struc 
ture, we present a volume rendering plot in Figure 
showing the magnitude of the current function A as va 
ume opacity and two sets of flux tubes (r ov = 1 and 0.9) 
as color lines. This figure demonstrates the relation be- 
tween regions of strong current and magnetic field lines 
for a = 60° rotator. 

First of all, we see that LOFLs are associated with the 
strong current layer inside the LC. Note that the current 
layer has finite thickness, and it encloses the LOFLs. 
Some field lines with r ov < 1.0 and r ov > 1.0 also lie 
within the strong current layer. Moreover, due to mag- 



netic reconnection in the current sheet (outside the LC), 
field lines can be closed beyond the LC, and we observe 
that some outgoing field lines with r ov > 0.95 stall and 
turn inward beyond the LC. These field lines trace the 
current sheet outside the LC to a finite distance before 
turning back. This effect makes finding the true LOFLs a 
more complicated exercise than we perform in this paper. 

Secondly, we see that field lines with r ov < 0.95 are 
generally open in the physical sense. Red curves in Fig- 
ure[5]show some representative field lines with r ov = 0.90. 
These field lines lie at the edge of the strong current layer 
inside the LC, and follow closely the current sheet out- 
side the LC. The current sheet outside the LC separates 
field lines originating from different poles. 

We conclude that r ov is a good tracer of current in the 
FF magnetosphere. For a = 60° rotator, strong current 
is associated with 0.95 < r ov < 1.05, and field lines with 
r ov < 0.95 lie outside the current sheet. For other incli- 
nation angles, the critical value of r ov — 0.95 varies, but 
not significantly. 

3. EMISSION CALCULATION 

The FF magnetosphere has no intrinsic acceleration 
and emission since the condition E ■ B = is satisfied 
everywhere by construction. Therefore, we will calculate 
the pulsar gamma-ray emission based solely on geometric 
grounds. Specifically, the calculation involves a prescrip- 
tion of the emission zone and a method to determine the 
intensity and direction of emission inside the emission 
zone. 

3.1. Emission Zone Geometry 

Here we summarize the emission zone geometry pre- 
scriptions by different theoretical models that we con- 
sider in this paper. 

The two-pole caustic (TPC) mo del is an exte nded ver- 
sion of the slot-gap model (e.g., Arons 1983). In this 
model, the emission is assumed to originate on LOFLs 
(i.e., r o y = 1) from the NS surface up to som e cutoff 
radius pyks fc Rudak|[2003| [Dyks et al.|[2004| . We use 
two cutoff radii: spherical radius r max , measured from 
the center of NS, and cylindrical radius R max , measured 
from the rotation axis. Together, these two parameters 
define the extent of the TPC emission zone. 

In this paper, we also propose a new emission model 
which we term the "separatrix layer" (SL) model. In this 
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Fig. 5. — 3D volume rendering of the current to flux ratio A, 
defined in equation |2jl, for the FF simulation with a = 60°. 
Transparency is determined by |A|, where only regions with large 
|A| are opaque. The color scale is the same as in Figure [II 
Blue curves indicate the LOFLs. while the red curves are field 
lines with r ov = 0.90. The two panels show different view- 
ing angles. Axes are in units of Rlc- The movie version of 
this plot is available in the online version of the paper and at 
http: //www. astro . princeton . edu/~anatoly/ cur sheet .mpeg 

model, we assume that emission comes from a layer in 
the vicinity of the separatrix, which separates open and 
closed field lines as well as f ield lines with different polar- 
alities. As we have seen in §2.3| the location of the sepa- 
ratrix corresponds to strong current layer/current sheet. 
This separatrix layer is described by r ov ~ 0.90 — 0.95, 
going from the stellar surface up to a cutoff radius in a 
way that is similar to the TPC model. The cutoff radii 
for SL model are not limited to be inside the light cylin- 
der. We will discuss the physical justification for this 
model in SjT] 



For the outer gap (OG) model, the emission is assumed 
to originate on the open field lines beyond the null charge 
surface (NCS), and extends up to the LC. Recent refine- 
ments of the OG models allow the inner bou ndary of 
the gamma-ray emission to be inside the N CS (|Hirotani| 



2007| [ Takata fc Chang|[2007| |Takata et aLl[2007[ p00«)- 

In order to incorporate such possibility, we assume that 
the emission comes from all field lines in a given flux 
tube that cross the NCS, extending from the stellar sur- 
face to certain cut-off radius. This approach maximizes 
the volume of the emission zone. The emission from the 
OG model is often assumed to be centered at certain r ov 
(e.g., r ov = 0.90). 

Unlike the OG model, the inner annular gap (IAG) 
model assumes that particles are accelerated in an an- 
nular open field li ne region extending fr om the NS sur- 
face to the NCS dQiao et al.||2004| |2007| ). The emission 
zone is bounded by LUFLs and the critical field lines, 
which cross the LC at the NCS. To simplify, we may also 
think of emission from the IAG to be centered at certain 
r ov < 1. 

For all models mentioned above, we consider the emis- 
sion originating from annular regions centered at fixed 
r nv w i th a G aussian width Ar ov = 0.025 similar to |Dyks| 
et al. (2004) as our standard emission zone prescription. 
The radial extent of the emission zone is controlled by 
the location of the NCS (for OG and IAG models), as 
well as the i? max parameter (we do not use r max in this 
paper) . 

3.2. Calculation Method 

In calculating the gamma-ray light curves we adopt 
two basic assumptions. First, we assume that the emis- 
sion is produced by outgoing energetic particles in the 
emission zone. These particles are confined to travel al- 
most along the magnetic field line, but may acquire some 
pitch angle 9 p . Such pitch angle is expected from both 
conventional SG and OG models, e.g., acqu ired during 



the pair creation process (Tang et al. 2008), or via the 



cyclotron resonant absorption (Harding et al. |2008| . If 
emission orig inates from the current sheet, as we postu- 
late in §7.2| we would expect even larger pitch angles 
since in the current sheet, magnetic field is relatively 
weak while the plasma is likely to be hot due to reconnec- 
tion. The emission mechanism can involve synchrotron, 
curvature, and inverse-Compton, but the direction of the 
emitted photon is along the direction of particle motion 
due to relativistic beaming. Calculation of the direction 
of emitted photon in the lab frame (LF), including the 
aberration and time delay effects, is discussed in detail 
in BS10. Here, we generalize these results to incorpo- 
rate the effect of finite particle pitch angle. Consider 
a drift frame (DF) in which the electric field vanishes. 
Let the velocity of the DF relative to the LF be Vd and 
Pd = Vd/c. In this frame, let 9 p be the particle pitch 
angle. Then the emission direction is determined by 

e = (a cos p )t+ (a sin 9 p cos <p)n+(a sm9 p sin Lp)b + (3d , 

where t, n, b are the tangential, normal and binormal 
unit vectors to the magnetic field line (in the LF), and a 
is determined by requiring that |e| = 1. 

This formula is es sentially the same as equation (6) of 
Takata et al. ( |2007[ ), except that it applies in any drift 



8 



Bai & Spitkovsky 



frame. The drift velocity can be chosen to be either the 
corotation velocity/?o , or the Ex B drift velocity. When 
9 p = 0, equation (nj) ensures that these two choices pro- 
duce the same result. However, this is no longer the case 
for finite pitch angle. This is because pitch angle is not 
Lorentz invariant and depends on frames. Moreover, the 

drift frame velocity with (3d = Po is no longer valid at 
R > Rlc since the frame velocity becomes superlumi- 
nal. Therefore, a better choice is to use the E x B drift 
velocity, which is the minimum drift velocity among all 
DFs, and this is the drift velocity we adopt throughout 
this paper. In our calculation, we further assume that 
the pitch angle is constant for all particles in the emis- 
sion zone. 

Our second assumption is constant emissivity along 
particle trajectories in the emission zone. In reality, as 
long as the particle emissivity does not vary dramati- 
cally in the emission zone, our calculation will be able 
to catch the overall features of gamma-ray light curves. 
Moreover, as pointed out in BS10, the overall appear- 
ance of the emission sky map is more sensitive to field 
structure, than to emissivity, because the formation of 
caustics largely depends on the details of field configura- 
tion. Adopting constant particle emissivity will help us 
locate possible emission zones in FF field by comparing 
with observations. For example, comparing model light 
curves with observations can help identify the radial ex- 
tent of the emission zone (i.e., -R max , see Sj5]), which then 
places constraint on the physics of radiation mechanism. 

One basic differ ence in our approach compared to pre- 
vious work (e.g., |Yadigaroglu||1997[ |Cheng et al.||2000[ 
Dyks et al.1|2004l 1'1'akata fc Chang||2007| |Harding et al | 



2008 ; Watters et al. 12009 ) is that the emission is weighted 



by the length along particle trajectories rather than along 
magnetic field lines. We've discussed in BS10 that parti- 
cle trajectory follows magnetic field lines in the corotat- 
ing frame (CF), and the magnetic field in the CF is the 
same as lab frame (LF) magnetic field. Therefore, to cal- 
culate particle trajectories, we trace magnetic field lines 
in the LF, and add the effect of stellar rotation. In Ap- 
pendix [O we provide a simple relation between segment 
length along a magnetic field line and along the corre- 
sponding particle trajectory. The difference is small well 
inside the LC, but the correction becomes important near 
the LC and beyond. Since FF field does not break down 
at the LC, we will allow the emission zone to extend 
beyond the LC, where such correction is necessarjQ 

Finally, in our calculation photons from the emission 
zone are projected to the sky map ((/>, £ bs), where 4> cor- 
responds to the rotation phase, and £ b s denotes the ob- 
server's viewing angle. A time delay correction is applied 
during the projection as usual 



e/R 



LC 



(6) 

The brightness of the sky map is proportional to the 
number of photons projected to the sky map bins, di- 
vided by the solid angle of each bin. The solid angle of 
each sky map bin is dfl = sm£ b s df[obsd<f>; therefore, the 
sky map brightness is proportional to photon counts in 

5 For example, consider a rotating split monopole field, where 
particle trajectories are straight lines, while the field lines are spi- 
rals. The FF field asymptotically approaches the split monopole 
field, and such correction avoids overweighing far-field emission. 



each bin divided by sin£ b s - Cutting the sky map at £ bs 
produces the light curve seen by the observer. The cor- 
rection of 1/ sin £ bs only affects the normalization of the 
light curves at each observer viewing angle, but not the 
shape of the light curves. Including it is important for 
the calculation of flux correction factor used to infer the 
intrinsic luminosity of gamma-ray pulsars. 

4. THE TWO-POLE CAUSTIC MODEL 

In this section we consider the emission from TPC 
model (r ov = 1.0) in FF field. The 3D shape of the 
last open field lines and the corresponding sky maps are 
shown in Figure [6] for inclination a = 60° . In order to 
better see the origin of the sky map features, we mark 
different field lines with different colors in panels (a) and 
(b), and use the same color convention to plot the sky 
map in panel (c). To guide the eye, we also plot six color 
rings crossing the field lines in panels (a), (b) and their 
projections onto the sky map in panel (c). Points on 
each ring are located at fixed distance, s, from the cen- 
ter of NS, measured along the field lines. The six rings 
span from s — 0.25i?LC to 1.50Rlc with an interval of 
0.25i?Lc- As the size of the star in the simulation is 
relatively large, R% m = 0.19Rlc, tracing emission from 
the surface out would leave large gaps on the sky map 
surrounding the poles. To fill these gaps, we add trac- 
ing of radiation from inside the star, where we impose a 
vacuum dipole field down to the radius Rn = O.OIRlc- 
This makes for effective size of the star for the purposes 
of sky map calculation to be R n- As FF field is nearly 
dipolar close to the star, the matching to the vacuum field 
is easy to achieve. In panel (d), we make the line plot 
of the sky map as in panel (c) , but the lines are colored 
with the value of 4-current length g, defined in equation 
Q. Since current drops rapidly with radius, the actual 
color scale is set by gr 2 . Note that in this panel, g = 
corresponds to the middle of the color bar. The black 
semi-circles around the pole indicate the size of the star 
in the simulation R% m . At smaller radii, where the FF 
field is not available, g is taken to be zero. Panels (e) 
and (f ) display the brightness of the sky map. Note that 
panel (e) shows the sky map brightness from the north 
pole only, as in all previous panels, while panel (f) is a 
complete sky map containing both poles. In calculating 
the sky map brightness, the synchrotron pitch angle 8 p 
is taken to be zero. 

First of all, the pattern of the sky map from the FF 
field, as seen in Figure [6j is very different from the sky 
map from TPC model using the vacuum magnetic field 
pyks & Rudak|[2003l [Dyks et al |[2004) see also Fig. M 
below). As noted in BS10, the aberration formula usecf 
in earlier works needs to be corrected, and this reduces 
the strength and increases the width of the caustics. The 
FF sky map here is also very different from the corrected 
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vacuum sky map (see Figure 4b, 4d of BS10, or Fig. 
below). The main differences include: 1) There is no 
caustic or local enhancement produced near the NS; 2) 
All caustics are formed near the LC. 

It is important to understand why the FF field gives 
such dramatically different appearance of the sky map 
even though the structure of the FF field is similar to 
that of the vacuum field close to the star. We find that 
this effect can be attributed to the shape of the polar cap. 
Recall that in Figure |3J the FF polar cap is more circular 
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Fig. 6. — Sky map and field line structure from the TPC model using FF field. The inclination angle is a = 60°. All field lines have 
r ov = 1.0 and are traced to the LC. a) and b): Field lines from the north pole viewed from two different directions, c): Projection of field 
lines from the north pole onto the sky map. Line colors are same as in a,b). Color rings represent points located at the same distance along 
the magnetic field lines from the center of NS: pink (0.25Rlc)> g re y (0.50i?Lc)i red (0.75Rlc)i green (I.ORlc), blue (1.25Rlc)> black 
(1.50Rlc)- Asterisks, circles, diamonds and stars are added to help in identifying field lines, d): Line plot of the sky map with each point 
colored by 4-current length q [see equation Q]; g = corresponds to the middle of the color bar. NCS is indicated by the grey dots, and 
the black circles correspond to the actual size of the star in the simulation (J?^ m ~ 0.19-Rlc). e): Sky map from field lines from the north 
magnetic pole, f) Sky map from both poles. In e) and f), zero brightness corresponds to the dark blue as indicated by "min" in the color 
bar. See text for more details. 



than its vacuum counterpart. More importantly, the FF 
polar cap is larger. The caustics in the vacuum field 
are significant in the open field line region, but weaken 
toward the LOFLs, tending to disappear if the emission 
is calculated from the field lines further in the closed 
zone. As the larger FF polar cap encloses the vacuum 
polar cap, the last open field lines of FF would map into 
the closed field zone of the vacuum magnetosphere. This 
explains the lack of caustics close to the NS surface using 
the FF field. 



Next we discuss the structure of the FF sky map for 
the TPC model (Fig. [fjj. Each pole contributes to two 
main caustic structures. For the field lines coming from 
the north pole (Fig. |6fc-e), the location of the caustics 
roughly coincides witn the green ring, and in the sky 
map it appears as a broad arc (Fig. |6^). At early phases 
{4> < 120°), the caustic is mainly caused by the sweep- 
back of field lines on the sky map (yellow field lines in 
Fig. [6^-c). This caustic is relatively weak and narrow. 
At later phases, there is a bright "knot" where emis- 
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Fig. 7. — The light curves (right) predicted from the TPC model 
using the FF field at various magnetic inclination angles a and 
observer viewing angles f b s - The corresponding sky maps are 
shown on the left. Emission zone is centered on LOFLs, extending 
to LC. A synchrotron pitch angle sin0 p = 0.1 is applied. The 
curves at the same a have the same normalization. 

sions from neighboring field lines cross each other on the 
sky map (near the intersection of light blue field lines 
and the red or green rings in Fig. pp). This knot and 
the adjacent region contribute to the most of the second 
caustic, which is brighter and wider. The knot can also 
be seen in panel (a), where light blue field lines fall close 
to each other. There is large shear between neighboring 
field lines, corresponding to large current flow. The field 
lines resulting in the knot enter the current sheet, and 
this feature persists in many different simulation^ 

Adding the emission from both poles, there canbe up 
to four peaks in the gamma-ray light curve. A large 
fraction of the area of the sky map includes contribu- 
tions from both poles. Besides the two bright caustics, 
emission from other regions is mostly smooth, except in 
regions close to the LC [e.g., the tips of red and cyan 
lines in panel (a)]. Comparing with Figure p3l we find 
that these regions are deeply embedded in the strong 
current layer, where the field structure may not be well 
resolved in the simulation. 

Figure J71 shows a gallery of representative light curves 
predictedfrom the TPC model in FF field for different 
inclination angles and observer viewing angles. In this 
plot we have chosen sin6* p = 0.1 to show the effect of 
non-zero pitch angle. We see that this smoothes the sky 
map substantially, as can be seen by comparing the a = 
60° sky map in this figure with Figure [6] Nevertheless, 
the light curves generally exhibit irregular shapes with 
up to 4 peaks. The central part of the light curves is 
typically brighter because both poles contribute. The 



6 We note that Contopoulos & Kalapotharakos (2009) observed 
similar features as tney evolve tne system tor longer time. We've 
also checked simulations with still high er resolution, and with finite 
conductivity using the prescription by |Gruz inov (2007bj), and find 
that the "knot" appears in all simulations. 



caustic structure due to field line sweepback is smoothed 
and mixed with the overlapping region from another pole, 
while the caustic from the vicinity of the "knot" is still 
prominent. Also note that some of the troughs in the 
light curves are due to the finite size of the poles. 

Comparing the light curves of Figure [7] with observa- 
tions, one generally fails to obtain a clean double-peak 
profile. At larger inclination angles, the trend is to have 
two very broad humps instead of narrow peaks. We have 
tried to restrict the geometry of the emission zone by 
tuning the parameters r max and i? maX i and found that 
they do not improve this situation. Moreover, we have 
also considered simulations with higher resolution and 
different diffusivity; the appearance of the sky map and 
the bulk behavior of the light curves does not change 
qualitatively between different runs. 

5. THE SEPARATRIX LAYER MODEL 

The difficulties in obtaining double-peaked light curves 
from last open field lines with the two-pole caustic model 
in FF field caused us to consider different locations for 
the emission region. In this section, we study the emis- 
sion from a set of field lines in the open volume that 
are still close to the LOFLs. This is motivated by the 
distributed nature of the return current in oblique ro- 
tators. We will refer to this model as the "separatrix 
layer" (SL), as the emission region is concentrated in a 
layer in the vicinity of the separatrix. In practice, we 
choose r ov = 0.9 — 0.95 for the SL model. In Figure [8j 
we plot the field lines and the corresponding sky maps 
with r ov = 0.9, for inclination a = 60°. The emission 
zone extends from the NS surface to cylindrical cut-off 
radius i? m ax = l-2i?z,c- We organize the six panels in 
the same way as in Figure [6j 

The SL model produces two bright and narrow caus- 
tics on the sky map. Comparing the TPC sky map with 
that of the SL model (i.e., moving from r ov = 1.0 to 
r ov = 0.9), we see that the weak caustic due to field 
line sweepback on the left of Figure [6b no longer exists, 
while the caustic associated with the "knot" has evolved 
to a strong caustic spanning a large portion of observer's 
viewing angle. The "knot" is no longer present. Com- 
bining the contribution from the two poles gives rise to 
two caustics in panel (f). 

The origin of the caustics can be traced from panels 
(a) to (c): we clearly see that the red, green, blue and 
black rings overlap at the locations of strong caustics. 
Therefore, the formation of the caustics is not due to the 
coincidence that emission from different field lines con- 
gregates on the sky map (as is the case for caustics in the 
vacuum field form), but the emission from one field line 
arrives simultaneously, piling up on the sky map. We 
call this effect "sky map stagnation" (SMS). In Figure[8j 
we see that SMS starts to develop at a distance from the 
star along field lines of 0.75Rlc (beyond the red ring). 
This corresponds to the outer magnetosphere. The SMS 
phenomenon covers about one half (ir radians) of mag- 
netic azimuth on each pole. Therefore, the full caustic 
structure caused by SMS sweeps a large fraction of the 
sky map, giving rise to two sharp peaks in the resultant 
light curve. 

In Appendix [D] we show that SMS is a natural conse- 
quence of the rotating split monopole solution of Michel 
( |1973b[ ). The direction of particle motion in this field 
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Fig. 8. — Same as Figure [6] but for the SL model, with inclination angle 
cylindrical radius R = 1.2i?£ q. 



60°. Field lines all have r Q 



0.9 and are traced to the 



is exactly radially outward, and the winding of the field 
lines into spirals compensates for th e time delay effect. 



The F F field of the aligned rotator (Contopoulos et al. 
1999 1 approaches the split monopole solution beyond the 
L(J. The oblique rotator also asymptotes to the inclined 
split monopole solution ( Bogovalov|[l9 99 , S06), and this 
causes the caustics in a general FT field. To check this, 
we have examined the field geometry of open field lines 
near and beyond LC. We find that the field resembles 
well the split monopole field. The deviations, which 
are stronger for lines with larger r ov , prevent SMS from 
occurring on all open field lines. We note, that split 
monopole field is just one out of many field line con- 
figurations that can result in sky map stagnation (see 
Appendix |D| . 



We also constructed sky maps from different open vol- 
ume coordinates, ranging from r ov = 0.95 to r ov = 0.60. 
We find that for smaller values of r ov the SMS behavior 
is stronger. At r ov = 0.60, SMS is achieved for all field 
lines. This suggests that as we shift towards the center of 
the open flux tube, the split monopole solution provides 
a bett er a pproximation of the field geometry (see also 
Figure [To| . We also examined the sky map from other 
inclination angles ranging from a — 15° to 90° (see also 
Figure [9] below). We find that SMS is in fact a general 
feature of the FF field. 

One consequence of SMS is that the intensity of the 
caustics is proportional to the length of particle trajec- 
tories within the emission zone, if emissivity is constant 
along the trajectory. We have chosen to cut off the ra- 
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diation at i? ma x = 1.2Rlc, which is arbitrary. Choosing 
larger cutoff radius will result in brighter peaks, and vice 
versa. As emphasized in ^j3j we have used the assumption 
of constant emissivity along particle trajectories rather 
than field lines. The two different treatments do not pro- 
duce significant differences at radii of order of Rlc ( our 
approach gives slightly fainter peaks); nevertheless, this 
approach is more self-consistent, and sets the framework 
for future work (e.g., spectral calculations). 

Another consequence of SMS is that the curvature ra- 
dius i? cur of the particle trajectories becomes much larger 
than the curvature radius of magnetic field lines. One 
extreme example is the curvature radius of the trajec- 
tory in split monopolar field, which equals infinity, even 
though the field lines are curved. We have calculated 
i? cur along the open field lines of the FF field (see Ap- 
pendix [C] for methodology) , and found that i? cur gener- 
ally increases as particle travels farther along the field 
line, and can reach about IORlc at the LC. Such large 
curvature radius makes curvature radiation potentially 
inefficient, and a non-zero pitch angle helps with radia- 
tion efficiency, which justifies our inclusion of this effect. 

We commen t on the inn er annular gap (IAG) model 

1 



by Qiao et al. 



fl2004 | 

n the 



2007| ) . The choice of r ov 
AG region 



0.90 is 

However, the emission 



consistent wit 

of the IAG model is thought to originate relatively close 
to the NS surface, up to the null charge surface (NCS). 
As seen from Figure [6j no caustic can form in the inner 
region of the magnetosphere since the field line is still far 
from approaching the asymptotic split monopole regime. 
The IAG model also suffers from the more serious prob- 
lem that only a small fraction of field lines cross the NCS 
(see fjg). 

Figure [9] collects the light curves of the SL model for 
five different inclination angles and four observer view- 
ing angles. The emission is centered on r ov = 0.9 (for 
a = 30° and 45°) and at r ov = 0.95 (for a = 60°, 75° 
and 90°). The emission is prescribed up to cylindrical 
radius R max = 1.5Rlc- Due to SMS, the intensity of 
the calculated light curves depends on how far the emis- 
sion is traced. To examine the relative importance of 
emission contributed from different regions, we also plot 
the light curves originating from within cylindrical ra- 
dius -R m ax = 0.9i?Ld indicated as grey lines. There are 
several features worth noting as we discuss below. 

First, we see that the double peak structure is a com- 
mon feature for most of the light curves. Similar to the 
conventional TPC and OG models using the vacuum field 



(Watters et al. 2009), the double peak feature is most 



likely to appear at relatively large a and Cobs- Second, 
we see that emission outside the LC dominates the peak 
brightness, while the bridge emission is from inside the 
LC. This is expected from SMS. As a consequence, the 
intensity contrast at peak and bridge reflects how far the 
separatrix layer extends beyond the LC. For many view- 
ing angles, the first peak tends to be stronger than the 
second peak. This is due to our assumption of constant 
emissivity and simplified choice of emission zone bound- 
aries. More physics has to be included to accurately con- 
strain the peak intensities. Finally, without smoothing 
due to the non-zero pitch angle, the SMS tends to pro- 
duce reflection-asymmetric light curves. That is, both 
peaks have a sharp rise and a relatively smooth decay. 
This appears to be inconsistent with some of the ob- 



served gamma-r ay pulsar light curv es [e.g., light curve of 
the Vela pulsar ( Abdo et al.||2009b ), which has a "horn" 
structure, symmetric upon reflection around the mid- 
dle of the bridge]. By adding a small pitch angle, the 
force-free light curves are smoothed and become more 
reflection-symmetric (e.g., the light curve with a = 60° 
and £ obs = 50°). 

It remains to explain why we choose r ov between 0.90 
and 0.95, justifying the emissi on zone as being in the 
"separatrix layer" . In Figure [Toj we show a series of 
sky maps with inclination a = 75° from r ov — 0.80 to 
r ov = 1.0, using sin# p = 0.1 smoothing. From the top 
panel, we see that at small r ov SMS develops on all field 
lines, but only covers a relatively small fraction of the 
sky map. Emission from two poles does not overlap. 
Consequently, 4 peaks appear in typical light curves. As 
one increases r ov , on the one hand, the sky map cov- 
erage from each pole increases, reducing the gap region 
between the two poles; on the other hand, the number 
of field lines that exhibit SMS decreases, reducing the 
intensity of two of the four peaks. At around r ov = 0.90 
to 0.95, emission from both poles starts to overlap on the 
sky map, and SMS remains strong for a number of field 
lines. Upon increasing r ov further towards last open field 
lines, a significant fraction of the sky map is contributed 
to by two poles, SMS weakens, and the prominent "knot" 
appears. A similar trend is observed for simulations with 
other inclination angles. 

In sum, we have shown that the separatrix layer model 
with r ov between 0.90 and 0.95 is capable of producing 
double peak light curves in a wide range of geometric 
parameters. We will discuss the origin of SL emission 
in fjTj The fact that the emission zone responsible for 
the caustics in the separatrix layer scenario is located in 
the outer magnetosphere is reminiscent of the outer gap 
model. However, there are two key differences between 
the SL model and the conventional OG. First, the OG 
model is based on the charge-separated magnetosphere, 
where the formation of the gap is thought to occur only 
on field lines that cross the NCS. On the other hand, the 
SL model assumes that the emission zone extends all the 
way from the stellar surface to beyond the light cylinder 
on all field lines at a given open volume flux surface, 
regardless of whether these lines cross the NCS. Second, 
at a fixed line of sight, most emission, especially the two 
peaks, are contributed by one pole in the conventional 
OG model. In contrast, in the above SL model, both 
poles contribute to the observed emission, particularly, 
with each pole contributing one peak. This is similar to 
the TPC model. 

6. THE OUTER-GAP MODEL 

In this section we discuss the sky maps and light curves 
from the OG model using the FF field. A key ingredient 
of the OG model is the null charge surface (NCS), which 
determines both the radiating field lines and the extent 
of emission on these field lines. In the case of vacuum 
dipole field, NCS is simply the surface where B z = 0. As 
a result, emission from the north pole mostly contributes 
to the 90° < Cobs < 180° portion of the sky map, and the 
south pole mainly contributes to the other half. In the 
FF magnetosphere, we can determine the NCS by finding 
where 47rp = V • E = 0. Guided by the recent work that 
proposed shifting the inner boundary of the OG towards 
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Fig. 9. — Same as Figure[7]but for the SL model. For inclination angles 30° and 45°, we've chosen r ov = 0.9, and the rest use r ov = 0.95. 
For all angles, sin# p = 0.1. tbr the sky maps and the light curves in black, the emission is traced to cylindrical radius R = 1.5-Rx,c- Grey 
curves show the contribution from the inner magnetosphere, R < OMR^c- 



the star (e.g., Hirotani|2007 ), we allow emission from the 
full length of each held line that crosses the NCS. This 
does not result in extraneous peaks as can happen in the 
vacuum field, but just maximizes the sky map coverage 
of the OG mod el, improving the bridge emission. 

In Figure |TT| we show the sky map of the OG model 
for FF field with r ov = 0.90. Because the separatrix layer 
model includes emission from all field lines with this r ov , 
the sky map of the OG model is simply a subset of that 
of the SL model. In fact, Figure[TT]is just part of Figure 
[9j except the value of r ov there is slightly different (0.9 
vs 0.95). Assigning emission only to NCS-crossing field 
lines almost always eliminates one peak from the light 
curve. The remaining peak is from the SMS effect of the 



south pole. The north pole contributes very little to the 
light curves. 

The reason for the elimination of one peak in the OG 
model is best seen in Figure [8] of the previous section. In 
panel (d) , the location of the NCS on the sky map from 
the north pole is indicated by grey dots. We can see that 
only a fraction of field lines cross the NCS. These field 
lines span roughly 180° in magnetic azimuth on the po- 
lar cap. This is to be compared with case of the vacuum 
field, as in Fig. |12[ where for similar geometry, almost 
all field lines cross the fi ■ B = surface (presumably 
the NCS) eventually. For the same r ov , the FF field lines 
are "straighter" than their vacuum counterparts, and a 
large fraction of them do not turn horizontal inside the 
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Fig. 10. — A comparison of sky maps with different r ov . The 
inclination angle is a = 75°. Field lines are traced to R = 1.2Rlc 



except for r 



1.0 which is traced to R = I.ORlc- sin#„ 



0.1 



is used for all plos. On the right shows the light curves with £ j-, s 
70°. 



e.g ., compare lines marked with circles in Figs 
12|. Therefore, using the FF field, half of the fiet 



LC 
and 

lines end up being non-emitting under the framework of 
the OG modeF] Consequently, half of the bright caus- 
tics in panel (ej of Figure [8] are missing. Combining the 
two poles, there is only one bright peak left. This is a 
dramatic difference from the results using the vacuum 
field. 



In Figure 11 we confirm that for all inclination angles, 



7 Note that the charge density inside the artificially large stellar 
radius of the FF simulation (-R|Jr m ) in Figure|s]is simply taken from 
the first term on the right hand side of equation n3B ; therefore, 
the grey dots inside the black circle are roughly a straight line at 
?obs = 90° (corresponding to B z = 0). Outside the black circle, 
the grey dots are generally below £ D bs = 90°, indicating that the 
location of the NCS in the FF field is further out than determined 
from the classical CI ■ B = criterion. The difference is more 
significant in regions marked with diamonds than regions marked 
with asterisks. Although there is a ju mp in the location of the 
NCS on the left of panel (d) in Figure [8] due to the large size of 
the star in the simulation, it is also clear irom the plot that it does 
not affect our selection of NCS-crossing field lines. 
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Fig. 11. — Same as 
zone is centered on r, 



Figure|9] 

3V = 0.90. 



but for the OG model, and emission 



most of the light curves have only one prominent peak. 
In many of these light curves, the shape of the peak as 
well as the bridge emission is very similar to those pro- 
duced by the OG model using the vacuum field, but the 
second peak is missing, or degrades into a weak hump. 
The weak hump is produced from the same pole as the 
main caustic, and mostly due to the continuation of the 
main caustic (e.g., the panel with a = 60°, £ h s = 50°). 
Given the weakness of the hump, it is unlikely to play an 
important role when more physics is added to relax the 
assumption of constant emissivity. 

In addition, we have experimented with the effect of 
inward emission. This e ffect was used in the original 
version of the OG model (Cheng et al.||1986a|b[ ), and re- 
considered in |Takata et al.| ( |2UU8| ). We find that the in- 
ward emission itselt will produce two peaks whose shapes 
are similar to the conventional OG model using the vac- 
uum field. Both peaks are relatively weak, and the bridge 
emission is almost half the brightness of the peaks. More 
problematically, the first peak lags the radio peak by 
more than 180°. Adding to this the main peak from 
the outward emission, there can be three peaks in total. 
Therefore, inward emission appears unlikely to play an 
important role. 

7. DISCUSSION 

We have shown that using the FF field, conventional 
TPC and OG model no longer work as well as in previous 
calculations using the vacuum dipole field. Alternatively, 
we find that the separatrix layer model, where emission 
comes from open field lines that lie just in the vicin- 
ity of the separatrix (including the strong current layer 
inside the LC and the current sheet outside the LC), is 
very promising in reproducing the prevailing double peak 
features of gamma-ray pulsar light curves. This model 
works due to the sky map stagnation effect, that is unique 
to FF field and occurs as the field becomes increasingly 
split-monopolar with radius. Our calculations of the SL 
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Fig. 12. — Sky map and field line structure using vacuum field presented for comparison with the SL model of the FF field from Fig. [9] 
The inclination angle is a = 60°. Field lines all have r ov = 0.9 and are traced to the LC (cylindrical radius < Rlc)- a ) an d b): Field lines 
from the north pole viewed from two different directions, c): Line plot of the sky map based on field lines from the north pole. Line colors 
are same as in a,b). Color rings represent points located at the same distance along the magnetic field lines from the center of NS: pink 
(0.25-Rto), grey (0.50Rlc)> re d green (1.0-R^c), blue (1.25-Rlc), black (1.50i?£c). Asterisks, circles, diamonds and stars are 

added to help in identifying field lines. Gray diamonds at (^ ^ a ~ 90° show the location of the null charge surface on these field lines, d): 
Sky map from the north pole that results from c) after including emission on every field line at r ov = 0.9. The emission is included from 
the star to the cutoff radius, e): Sky map that results when the emission from both poles is superimposed. This map incorporates emission 
regions of both TPC and OG models (formally, TPC model is usually centered on r ov = 1, but the picture is similar at r ov = 0.9). While it 
generally results in too many peaks, this plot demonstrates the origin of all caustics in the vacuum field. The clustering of rings in panel c) 
indicates enhanced caustic emission in the sky map, although here the caustics are caused by near overlap of neighboring field lines, rather 
than SMS from each field line, f): Sky map that results from e) when emission below null charge surface is removed. This is the classical 
OG model. 



model in Sj5]work mainly on the geometric grounds with 
little physical input. In this section we discuss the im- 
plications of these results to pulsar gamma-ray radiation 
mechanism, and compare our results with observational 
data. 

7.1. Applicability of the FF field 

The FF field provides a more realistic model for the 
pulsar magnetosphere than vacuum held because it takes 
into account the field distortions due to self-consistent 
currents caused by pulsar spin down. The basic assump- 
tion of the force-free approximation is that pair creation 
is so efficient that it fills the magnetosphere with abun- 



dant pair plasma to make it a perfect conductor every- 
where. Real pulsar magnetospheres are more compli- 
cated due to dissipation and reconnection in the current 
sheet, for mation of gaps, and the presence of differential 
rotation ( jTim okhin 2 007a|b I. 

One way of testing the robustness of the result ob- 
tained in this paper is to consider field s from strong-field 



electrody namics (SFE) introduced by Gruzinov (2007a 
|2008a|b| , which generalizes the FF equations to include 
dissipation. In SFE, the FF condition E B — is 
relaxed in space-like regions, where the Ohmic dissi- 
pation is present, while time- like regions are still non- 
dissipative. The dissipation and reconnection in the FF 
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current sheet may be better modeled with SFE. In addi- 
tion, vacuum gaps in pulsar magnetosphere may also be 
considered as dissipative regions, where E ■ B ^ 0. We 
have implemented SFE formulation into our FF code by 
replacing t he FF current by the SFE current [see equa- 
tion (6) of Gruzinov (2007b)]. The conductivity scalar 
a = <j(Eq, no) is taken to be constant. The natural scale 
of a is fl/c 2 , which we denote by a = 1. We have run 
our simulation with different values of a around 1. We 
find that the overall field configuration is more smooth 
than the FF field configuration. Particularly, the "knot" 
in the sky map from LOFLs which penetrate the current 
sheet becomes less prominent. The strong current layer 
inside the LC is not changed qualitatively, but the cur- 
rent sheet outside the LC is smoothed appreciably. More 
energy is dissipated in the magnetosphere, and the total 
Poynting flux decreases by about 10% over about IRlc- 
Despite the fact that SFE produces smoother magne- 
tospheric structure, we find that the resulting sky maps 
and light curves do not differ qualitatively from the re- 
sults with the FF field presented in this paper. There- 
fore, we believe that our result is robust and may also be 
applicable to dissipative magnetospheres with moderate 
levels of dissipation. 

7.2. The Origin of the Separatrix Layer Emission 

We have shown that if pulsar high-energy emission 
originates from open field line regions in the vicinity of 
the separatrix, extending from the NS surface to beyond 
the LC, two bright caustics will form due to sky map 
stagnation. Currently, there is no theory directly point- 
ing to the existence of the SL emission; however, based 
on the geometric location of the SL, there are clues about 
its origin. 

1. Modified SG/TPC model. We note that SMS is 
already significant at r ov = 0.95, which is very close to 
LOFLs. Also, in our calculation, the emission zone is 
assumed to have a Gaussian profile with width Ar ov = 
0.025 centered on r ov . Therefore, the edge of the emis- 
sion zone can be considered on the LOFLs. This emission 
zone geometry is consistent with the SG model. There- 
fore, if the thickness of the pair formation front is as large 
as Ar ov = 0.025 — 0.05, a modified SG model may be able 
to produce the desired two peak light curves. Such SG 
must extend at least to the LC to produce sharp peaks. 
Recent SG/TPC models consider em ission from smaller 
r ov flux tubes (e.g., |Venter et al.|2009 ). Given that in FF 
field only the field lines with r ov < 0.95 are truly open 
and do not close through the current sheet, it is likely 
that at least geometrically the SL model can be consid- 
ered as an extension of TPC and invokes the same field 
lines. The physical reason for the acceleration on these 
lines is unlikely to be a slot gap, however, and is more 
related to reconnection in the current sheet, as discussed 
below. 

2. Modified OG model. Another clue comes from 
the distribution of space-like regions in the magneto- 
sphere. Regions of space-like 4-current, which are likely 
to develop instabi lities and could have accelerating fields 

Gruzinovj|2007b[ ) , are preferentially located beyond the 
US. Similar locations are typically invoked in OG mod- 
els, although for reasons that rely on complete charge 
separation in the magnetosphere. Even if OG-type model 
can be justified based on the sign of 4-current, it would 



not be sufficient to produce two peaks - the emission 
has to come from more field lines than those crossing 
the NCS. Interestingly, looking at Figure [8]i, the regions 
showing the pile-up due to SMS have values of g near 
zero, or negative. If small positive values of g can also 
cause acceleration, this could explain the emission from 
the whole flux tube. 

3. Association with the strong current. The presence 
of a strong current layer within the LC and a current 
sheet beyond LC are general features of FF field. Strong 
dissipation and reconnection is exp ected in such regions 
(Gruzinov 2007a Lyubarsky 2008), and the energy re- 
lease can reacfi a tew percent of the spin d own power, 
proba bly radiated in X-rays and gamma rays ( |Lyubarskii 
19961. Using Bogovalov's (1999 ) current sheet so lution 
tor trie inclined split-monopole, [Kirk et aT| ([2002 ) were 
able to reproduce the typical double-peak light curves, 
indicating that the field configuration around the current 
sheet beyond the L C is also geometrical ly fa vored. Based 
on the same model, |Petri fc Kirk| ( [2005] ) and |Petri|(|2009 | 
further reasonably well reproduced the optical polariza- 
tion of the Crab pulsar and phase-resolved spectrum for 
the Geminga pulsar. Our result that the SL emission 
zone resides in the vicinity of the strong current layer 
and/or current sheet (see Figure [5| implies a connection 
with the strong current. At higher inclinations, the re- 
turn current becomes increasingly distributed over the 
open flux, and it is possible that some of the accelera- 
tion takes place on the open field lines. Also, the caustics 
of emission form in the outer magnetosphere, and the ac- 
celeration on these field lines can be affected by the prox- 
imity of the Y-point. Reconnection near the Y-point (or 
"Y-ring" in 3D) can load the open field lines with plasma, 
and cause field -aligned acce lerating fields, akin to auroral 
double-layers ( | Arons|[2008[ ) . Emission from the equato- 
rial current sheet beyond the LC can also be important, 
as the field lines we invoke for SL trace close to the cur- 
rent sheet in the wind zone. Modeling emission from the 
current sheet itself is more complicated in FF simulation, 
since ideal MHD doesn't apply in the current sheet, and 
the magnetic field direction abruptly changes in the unre- 
solved current sheet, which leads to tracing errors. This 
will be studied more in future work. The facts that the 
bulk of the SL emission is formed at and beyond the LC, 
and that emission from two poles overlaps on the sky 
map and forms two peaks only close to the edge of the 
open flux tube (for r ov > 0.9, Figure 10), suggests that 



the current sheet beyond the LC is responsible for the 
large part of gamma-ray emission. 

Our discussion on the origin of the SL model is highly 
speculative. A detailed theory of the electrodynamics of 
current sheets has to be constructed to test these scenar- 
ios and uncover the origin of the SL emission. 

7.3. Comparison with Observations 

One year after launch, Fermi LAT has greatly ex- 
panded the number of detected gamma-ray pulsars either 
by identifying previously unresolved/unidentified sources 
(e.g., known radio pulsars, millisecond pulsars (MSP), 
EGRET sources, sup ernova remnants (SNR) pu lsar 



Abdo et al.|[2008l | 2009a| c]g|), or 
AT sources lAbdo et al.l|2009f). 



wind nebulae (PWN), 
by the blind search of 
Together with previously known gamma-ray pulsar data, 
a sample of more than 40 gamma-ray pulsars is now avail- 
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able ( Abdo et al.|2009"h ), allowing more systematic com- 
parison between observations and theoretical models. 

Under the framework of the SL model with FF field, 
the double-peak light curves can be produced over a large 
range of geometric parameters, as shown in ^5} As a ge- 
ometric model, it provides robust prediction on the tem- 
poral location of the two peaks regardless of the input 
physics. The phase lag of the first peak relative to the 
radio peak (which likely coincides with the location of the 
polar cap on the sky map) ranges from about 8 = 0.14 
to larger than 5 = 0.3, and the separation between the 
two peaks ranges from zero (i.e., one peak) to S — 0.5. 
Smaller phase lag typically corresponds to larger separa- 
tion, and vice versa, as seen from the caustic structure 
from Figure [9[ Also, larger inclination angles tend to 
result in larger peak separation. Statistically, our model 
predicts that the majority of light curves should have 
widely-separated double peaks, pulsars with smaller peak 
separation are less common, pulsars with only one peak 
are possible but rare. Another consequence of the SL 
model is that smaller peak separation usually indicates 
stronger bridge emission between the peaks. The best 
example of this can be drawn from the light curve with 
a-45 o ,U s = 50 o in Fig. § 

Based on the current data sample of gamma-ray pul- 
sars, the majority of them show widely separated double- 
peak pulse profiles, in good agreement with the SL 
model. A small fraction shows double peaks with smalle r 
separation [e.g., B1716-44 ( jPellizzoni et al.| |2009a|) , 
J0007+7303, J1459-60, J1741-0254 ( [Abdo et al'|[200<jf| )J, 
but with relatively strong bridge emission, also m qual- 
itative agreement with our result R A few objects show 
only one peak [e.g., J0357+32 dAbdo et afl [20091 , 
PSR 0437-4 715, PSR 0613-0200 Abdoet^ [20094, 
J2229+6114 flPellizzoni et al.|2009b Ij. Such pulse profiles 
may be obtained when the observer s line of sight tangen- 
tially cuts through the caustics on the sky map (examples 
include a = 30°,£ o b s = 50° in Figure [9]). The small num- 
ber of such pulsars agrees with the small chance of such 
tangential incidence. Noticeably, the fraction of single 
peak MSPs appears to be highe r than for normal gamma- 
ray pulsars ( |Abdo et al.|2009g ) . This may reflect the fact 
that MSPs tend to have smaller inclination angle which 
increases the chance of tangential incidence. Given the 
relatively old age of the MSP population, their smaller 
inclination angle might be caused by alignment torque s 
pavis fc Goldstein|1970[|Goldreich|1970 [MelatosT 2000). 



The distribution ot phase lags from the SL model is 
not directly comparable with the latest data since most 
new gamma-ray pulsars are radio-quiet. For available 
radio-loud pulsars, we find some agreement, but not al- 
ways. For example, the phase separation between peaks 
in B1706-44 appears to be too small to fit our model, as 
both peaks come in the first half of the rotation. There 
may be an additional phase lag from the radio peak due 
to high altitude radio emissio n from energetic pulsars 
(Weltevrede & Johnston 2008). Therefore, it is uncer- 
tain whether observed phase lag can be used as a reliable 
diagnostics. 



8 Note that in Figure 2 of | Abdo et"al] ( |2009f| ) , the first peak is 
placed at phase of 0.3 because most ot the new gamma-ray pulsars 
are radio quiet. 



8. SUMMARY AND CONCLUSION 

This paper is concerned with high-energy emission 
from gamma-ray pulsars modeled by the force-free (FF) 
field. The FF field takes into account the effect of con- 
ducting plasma, and the charges and currents are solved 
self-consistently using time-dependent simulations. By 
construction, the FF magnetosphere does not emit be- 
cause E ■ B = everywhere, but it provides a more 
realistic approximation to the field structure. The FF 
field differs substantially from the vacuum dipole field in 
several ways. First, there is a larger magnetic flux on 
the open field lines. As a result, the polar cap is larger 
than the vacuum polar cap. Second, beyond the LC, FF 
field has a strong current sheet and field lines approach 
inclined rotating split-monopole shape. Inside the LC, 
there is a strong current layer coinciding with the last 
open field lines. This layer is a thin current sheet for the 
aligned rotator, which connects to the current sheet out- 
side the LC. The amount of current carried by the current 
sheet inside the LC decreases with the inclination angle. 
The inner current sheet disappears completely for the or- 
thogonal rotator. Current distribution on the polar cap 
varies with a, with larger fraction of the current closing 
through the polar cap at larger inclinations, rather than 
in a thin boundary layer. Third, the charge distribu- 
tion in the FF magnetosphere contains significant con- 
centration of charge in the current sheet /strong current 
layer, and the charge distribution in the current sheet for 
the aligned rotator is substantially different from that of 
the oblique rotators. Moreover, the null charge surface 
(NCS) is no longer determined by ft ■ B = 0, and a large 
fraction of open field lines does not cross the NCS at 
all. Fourth, the 4-current in substantial fraction of the 
FF magnetosphere is space-like, requiring the presence 
of charges of both signs, which is inconsistent with com- 
pletely charge-separated picture. 

The differences in the field structure of FF and vac- 
uum fields give rise to substantial differences in resulting 
theoretical light curves. The polar cap in the FF mag- 
netosphere is more circular and larger than the vacuum 
polar cap. Therefore, close to the star, the location of 
LOFLs of the FF field corresponds to closed field line 
regions in the vacuum field. Due to this effect, the caus- 
tics seen in the two-pole caustic (TPC) model using the 
vacuum field no longer exist in force-free field. Instead, 
the pattern on the sky map is more complicated: up to 
four peaks (caustics) may exist, and a large fraction of 
the sky map is covered by emission from both poles. In 
addition, there is a bright "knot" region on the sky map 
associated with the strong current. In any case, we find 
it difficult to reproduce the typical double-peak profile 
in the resulting light curves using the TPC model. 

We have shown that if we relax the emission zone in 
the TPC model into a separatrix layer (SL), where the 
emission comes from the open field line regions near the 
strong current layer/current sheet and extends from the 
NS surface up to and beyond the LC, two bright caus- 
tics will appear. The cause of the caustics is completely 
different from that in the case of a vacuum field. For 
a large fraction of open field lines in the outer magne- 
tosphere, the emission from different heights along the 
field line piles up at the same spot on the sky map, which 
strongly enhances the sky map brightness at that loca- 
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tion. This effect, which we call "sky map stagnation" 
(SMS), is due to the fact that the FF field asymptoti- 
cally approaches the split monopole field at large radii. 
As a result, the intensity of the caustics due to SMS effect 
is determined by the extent of the emission zone. The 
typical feature of the light curves from this SL model 
is two sharp peaks with more or less uniform off-peak 
emission. Each pole contributes to one peak. The rise of 
the peak is usually sharp while the drop is more smooth. 
Adding the effect of non-zero pitch angle smoothes the 
sky map. This model provides robust prediction on the 
location and the separation of the two peaks, and the 
resulting light curves match the general features of most 
observed gamma-ray pulsar light curves. Although there 
is no theory predicting the existence of this separatrix 
layer emission, its location in the pulsar magnetosphere 
strongly suggests its association with the strong current 
layer and/or current sheet, particularly the current sheet 
beyond the light cylinder. 

The sky map from the OG model can be considered as 
part of the SL model, restricted only to field lines that 
cross the NCS. We find, however, almost always that 
only one peak can form. The reason is that a substantial 
fraction of the open field lines do not cross the NCS at all. 
The coverage of the emission on the sky map is similar 
to that of vacuum field, but the formation of caustics is 
due to SMS. Some gamma-ray pulsars recently observed 
by Fermi telescope show only one prominent peak, and 
the OG model using the FF field may be relevant to 
such pulsars. The "inner annular gap" (IAG) model also 
relies on field lines that cross the NCS, and thus suffers 
the same problem as the OG model. However, because 
the emission zone in the IAG model is relatively close to 
the NS, where SMS does not develop, it will be even more 



difficult to produce any peaks with this model using FF 
field. 

In constructing models of the emission zones, we as- 
signed emission to flux tubes concentric to the last open 
field lines. This can be an oversimplification, and the 
actual emitting flux tube can have a more complicated 
shape. This is particularly so if reconncction in the Y- 
point region affects neighboring field lines in the outer 
magnetosphere. The origin of those field lines on the po- 
lar cap does not have to be concentric to LOFLs. An 
avenue that we have yet to explore in detail is to as- 
sign emission not per field line, but in the whole volume 
of the magnetosphere, weighted by some proxy of the 
magnetospheric dissipation, such as the strength of the 
current. The gross features of the force-free sky maps, 
such as the caustic formation in the outer magnetosphere 
due to sky map stagnation will likely remain the same; 
however, quantitative shape of the light curves may be 
affected. This will be the subject of future work. In 
addition, our modeling of pulsar's high-energy emission 
assumes constant emissivity along particle trajectories, 
which works on purely geometric basis and serves as the 
first step in comparing different theoretical models and 
identifying possible regions of emission. More physics 
input is needed for a consistent study of pulsar's high- 
energy emission that would reveal the origin of pulsar 
magnetospheric accelerator. 
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APPENDIX 



A. DERIVATION OF USEFUL RELATIONS IN FORCE-FREE MAGNETOSPHERE 



Formal proof of equations ([I]) and (p| can be found in Gruzinov (2006). In this appendix we summarize Gruzinov's 
proof and further explain the physicaTmeaning of A. 

Both equations follow from two assumptions: a) The electromagnetic field is force-free; b) The electromagnetic 
pattern corotates at angular velocity £1. Condition b), when applied to an arbitrary vector field U , can be written as 



[d t + (n x r) ■ v][/ = n x U . 

For electric and magnetic fields, after some algebra, the above equation can be reduced to 



dE 

~dt 



V x [(SI x r) x E] + Airp(Sl x r) 



(Al) 



(A2) 



dB 

~dt 



= V x [(SI x r) x B] 



(A3) 



In line with 



Gruzinov 



where x is a scalar tunc 



(2006), equation (A3 1, together with the induction equation, implies E = — /3 x B + Vx, 
ion to be determined. The fact that the NS is a perfect conductor ensures x to be constant 
on the stellar surface, while the FF condition requires x to be constant along magnetic field lines. Therefore, x must 
be constant everyw here, completing the proof for equation (JT]). 

Gruzinov ( 2006 ) proved equation ^ by using the variational principle. Despite the mathematical beauty of this 
method, the physical meaning of A is not explicitly seen. Below we prove equation (pj) from direct calculation. 
We begin by substituting equation ( |A2[ ) into the Ampere's law of Maxwell's equation, and making use of equation 
1|), one immediately obtains 



(1) 



Vx[B + ftx(ftx B)} = j — p(Sl x r) 



(A4) 
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It re mains to show that B \\ (j — p(f2 x r)). This becomes obvious since the FF current can be cast into [cf. equation 
(2) of |Spitkovsky] poSg] )] 

Ex B 

J = P — 



B 2 



(A5) 



where jy is the current density that is parallel to B. Because of equation (|1J), the corotation velocity fl x r can be 
decomposed into the E x B drift velocity and a component that is parallel to B. 
In sum, we have reached a remarkably simple and clear formula 

XB=j p(n xr) = f , (A6) 

where j c is the current density in the corotating frame (CF) ( Schiff|l939 Gr0n|l984 ). Therefore, in CF, the magnetic 
field line and current are parallel to each other. The ratio between j u and B is denoted by A, and is constant along 
the magnetic field lines. 

B. CURRENT SHEET IN THE FORCE-FREE FIELD 
In this appendix we show that the cu rrent in the curren t sheet outside the LC in FF field does not necessarily connect 



to the star. Our discussion adopts the 
monopole, and asymptotically approac 



Bogovalov] ( 1999[ )'s solution. It describes the field structure of a rotating split 
les to the FF field at large radii. The Bogovalov's solution reads 



B, 



fo i 



Re^r, 6,<t>) , (Bl) 

where the subscript "to" represents monopole, fo is an arbitrary constant, R = rsin# is the cylindrical radius, and 

rj(r, 9, </>) = Z?(sinasin6>sin(</> + r/R^c) + cosacosfl) = D(x) , (B2) 

where D{x) = 1 for x > and D{x) = — 1 for x < 0. 
The current sheet lies on the surface determined by x = 0, namely, 



cot 9 = tana sin(</> + r/R^c) ■ 
Clearly, we have ir/2 — a < 9 < ir/2 + a. One can find the normal direction of this surface to be 

1 



n = 



(Ae r - eg + Ae 4 



where 



.4 



tan 2 a sin 2 ( 



r/Rhc)] tan a cos(^ + r/i? LC ) 



(B3) 



(B4) 



(B5) 



Th e natural frame to calculate the current in the current sheet is the CF, where the current is given by equation 
(A4|. In addition to the curl of B, the CF current has contribution from the curl of f3 x E term. Therefore, it is 
useful to define 



ft, 



B, 

Jo 
r 



n 

+ A) x (J3 x B rn ) 
2 [(l-i? 2 )e r -i?e ]77(r,0,, 



(B6) 



The surface current density in the current sheet in the CF is then given by 

J = n x (H+ - H~) , (B7) 

where superscripts "+" and "— " denote the quantities just above and beneath the current sheet. 

It is useful to look at two special cases, namely, the aligned and orthogonal rotators. For a — 0, we have n = e z . 
From the above equations, it is obvious that J has both radial and azimuthal components. The radial component 
connects the current to the star, and decreases as 1/r. The azimuthal component asymptotically approaches a constant. 
Integrating over a circle with radius r, one obtains the net current carried by the current sheet that connects to the 
star 



/„ = 



rd(j> = 47t/ 



(B8) 



Note that Iq is a constant that does not depend on r, reflecting the conservation of current. One can also check that lo 
equals to the amount of current flowing from the star (outside the current sheet), but with a minus sign, as expected: 



W,o = f (V x J?) 



e r r sin 



= -la 



(B9) 



Q=0 



For a — 90°, A —> oo, thus n = (e. r + e^)/y/2. Plugging into the above equations, one finds that J has only the 9 
component. This means that the current in the current sheet is not connected to the star, but forms poloidal loops in 
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Fig. 13. — The amount of current in the current sheet of the r otati ng split monopole that is connected to the star ishect as a function of 
mag netic inclination angle a. 7 s heet is calculated from equation (|BX0|> , normalized to the value for the aligned rotator Iq, given by equation 



the 6 direction. In the FF field from our simulations, we do observe such loops in the current sheet outside the LC. 
Since they are not connected to the star, it becomes clear that there is no current sheet in the inner magnetospherc 
for the orthogonal rotator. 

After some algebra, we obtain the general expression of the current carried by the current sheet that is connected 
to the star 

4heet (oQ = f 2 * sing dcj> 

lo ' ~ Jo VI + 2A 2 2tt ' 1 ' 



where 6 is determined by equation (B3| and A is given by equation (B5). The quantity does not depend on r, as 
expected. Moreover, one can check that this quantity also equals to the net current flowing from the the rest of the 
star 

/star ( a ) = — L I cos sin 0D{x)d0d(j) , (Bll) 



2tt 



where x is given by equation ( B2 ) . The dependence of / s hcct on inclination angle a is shown in Figure 



the amount of current in the current sheet that is connected to st ar m onotonically decreases with 
to the degradation of the current sheet inside the LC discussed in §2.1| 



rep] A 
a/Thi 



We see that 
his is related 



C. SEGMENT LENGTH AND CURVATURE RADIUS OF PARTICLE TRAJECTORIES 



Here we describe the method for calculating the length and curvature radius of the particle trajectories. For simplicity, 
we consider particles with zero pitch angle. The discussion here is valid so far as particle pitch angle is small, as we 
use in this paper. The direction of particle motion consists of a combination of E x B drift velocity (or corotation 
velocity) and velocity parallel to the magnetic field 



e = Pd + oct = A) + a't 



(CI) 



where t = ±B/B is the unit vector parallel to the magnetic field line (without loss of generality, hereafter we assume 
the radial component of t points outward), and a or a' are obtained by requiring |e| = 1: 



p d -t±d(p d -t)* + l-0>i , 



a 



/V*± vW*) 2 + l-/3 2 



(C2) 



where the plus/minus signs correspond to trajectories along/opposite to magnetic field line. Usually, we expect /3 -t < 
due to field sweepback. Therefore, the outgoing particles always have a' > 0, and are expected to have a' > 1 (for 

reference, split monopole field has a' = yl + (R/ Rlc) 2 )- For ingoing particles, however, a' is positive inside the LC, 
negative outside the LC. This means that ingoing particles can return to the star only if they are generated inside the 
LC. Below we consider only the outgoing particles, taking the plus sign in the equations above. 

Particle trajectories can be obtained by traci ng th e direction of motion while rotating the star at the same time. 
Since e has a corotation component, equations (CI) ensure that particle trajectories follow the field lines exactly in 
the corotating frame (CF). Although this is obvious, we still go through the proof, which will be useful for numerical 
calculations of curvature radius of particle trajectories. 

The proof is best demonstrated in a graphical view. In equation (CI ), the vector equation forms a triangle as shown 
in Figure [Mj The length of each edge is expressed in terms of velocity. Multiplying each of the edges by an amount of 
time dt, we obtain another interpretation of the graph. The upper edge gives an infinitesimal segment length dl = cdt 
along particle trajectory in the direction of e. The lower edge becomes ds = a' cdt, which marks the corresponding 
segment length along the magnetic field line. As the particle travels by dl, the star (and the field lines) rotate by 
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an angle d(f> = ildt — fldl/c = dl/R^c- Therefore, the entire field line is shifted by Rd(f> = P^cdt = Rdl/R^c in the 
azimuthal direction, this corresponds to the right edge. Now we've obtained a closed triangle that is similar to the 
original velocity triangle. Based on this new interpretation, we conclude that particle trajectory follows the magnetic 
field line exactly in the CF. 

The fact that particle trajectories follow CF magnetic field lines as well as the similarity relations from the velocity 
and segment triangles provide a natural way to reconstruct particle trajectory. The key relation is between dl and ds, 
which satisfies 

dl = ds/a' . (C3) 

In practice, we originally have a series of segment lengths along field lines s[i]. From this formula, one can obtain a 
series of segment lengths of the corresponding particle trajectory l[i]. This further provides the phase of rotation as 
the particle travels, cj>[i]. Finally, particle trajectory can be reconstructed by rotating each point on the magnetic field 
line by <j)[i]. It is then straightforward to (numerically) calculate the curvature radius of the particle tra ject ory. 

When performing the calculation, however, there can be cases where the square root in equation (C2| becomes 
imaginary. This always happens in the vacuum dipolc field once the particle travels out of the LC; it also h app ens 
occasionally in the current sheet since the ideal MHD doesn't apply there and one does not expect equation (CI) to 



hold. Moreover, across the FF current sheet, the toroidal component of the magnetic field changes sign, resulting in 
ambiguity in defining "ingoing" and "outgoing" directions. In su ch r egions, the proposed "outgoing" particles may 
have a' < 1, which is not desired. Therefore, we modify equation (C3) into 



dl = ds/ max(a , 1) 



(C4) 



In calculating a', if the terms under the square root is negative, we just take them to be zero. 

One can further show the analytical expression for the curvature radius of particle trajectory based on the discussions 
above: 



1 



fl x e 



, de 

ds 



(C5) 



Numerical experiments show that curvature radius obtained from this approach has very similar accuracy as the first 
method. 



D. CONDITIONS FOR SKY MAP STAGNATION 



According to the formula of aberration (CI I and time delay ([6]), the condition for emission along a magnetic field 
line to become stagnant on the sky map is 



B ■ Ve z = 0, 



B ■ V( 



e) = 



where 4> v denotes the azimuthal angle of the emission direction t). 
The second equation can be reduced to 



B ■ V(r • e) + [e X (B ■ V)e] 2 /(1 - e\) = . 



(Dl) 



(D2) 



Note that the emission direction depends only on the direction of the magnetic field, but not its strength. Therefore, 
the only variable left is the direction of the field, which has two degrees of freedom. The a bove are two independent 
equations for B. Therefore, the direction of B can be uniquely determined from equations (Dl), if initial condition is 
given. Nevertheless, the equations above are still too complicated to solve analytically. Instead of solving the equations 
directly, we test certain field geometrie s to see if the stagnation condition is satisfied. 



We consider the split monopole field (Bl ) and show that it does indeed satisfy equations (Dl ). According to equation 
(CI), the emission direction from split monopole is exactly along the radial direction e = e r . Since the strength 
fo/r 2 and sign r)(r,0,<f)) have nothing to do with the emission direction e m , we can safely omit it in the calculations 
below. 
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Now we substitute the split monopole field into these conditions. The fi rst e quation in (Dl ) becomes B 
Note that _B™ = 0; therefore, this condition is satisfied. For equation (D2), the first term is just B r ■■ 

,2 



"■Vcos(9 = 0. 
= 1, while the 
emission from 



second term can be shown to be equal to —Rsm8/r(l — sin^ 8) = —1, and they sum to zero. Therefore 
any one field line in the split monopole field is projected exactly to the same p oint on the sky map. 

Obviously, the split monopole is not the only possible solution to equations (Dl). It only corresponds to a special 
initial condition. In fact, we see in Figure [8] that at r ov = 0.9, SMS starts to develop from inside the LC, where the 
FF field still deviates from the split monopole solution. Nevertheless, SMS is much easier to achieve at smaller r ov 
(see Figure 10 ), or at larger radii (R Rlc)- The fact that the FF solution approaches the inclined split monopole 
asymptote (" Bogovalov||1999 ) , as well as the observation of SMS strongly suggest its connection to the split monopole 
field geometry at large radii. 
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